#include <iostream>
#include <fstream>
#include <cassert>
#include <string>
#include <vector>
#include <sstream>
#include <cmath>
#include <ctime>
#include <thread>
#include <mutex>
#include <future>
// 每个采样的行数
#define LINES 35305
// 文件头行数
#define H_LINE 9

using namespace std;
// 锁
mutex mu;
template <class T>
void time(T &&func)
{
    using namespace std::chrono;
    auto start = system_clock::now();
    func();
    duration<double> diff = system_clock::now() - start;
    cout << "Time:" << diff.count() << "seconds" << endl;
}
// 打印向量
template <typename T>
void show_info(vector<vector<T>> &v)
{
    int col = v.size();    // 行数
    int row = v[0].size(); //  列数
    for (int i = 0; i < col; i++)
    {
        for (int j = 0; j < row; j++)
        {
            cout << v[i][j] << " ";
        }
        cout << endl;
    }
}
// 计算距离
double dist(vector<double> &a, vector<double> &b)
{
    return sqrt(pow(a[0] - b[0], 2) + pow(a[1] - b[1], 2) + pow(a[2] - b[2], 2));
}
// 字符串分割解析
void split(const std::string &str, std::vector<double> &val, int count)
{
    auto pos = str.find(" ");
    if (pos == std::string::npos)
    {
        val[count - 2] = std::stod(str);
        return;
    }

    auto s = str.substr(0, pos);
    if (count > 1)
    {
        val[count - 2] = std::stod(s);
    }
    return split(str.substr(pos + 1), val, count + 1);
}
// 打印一维向量
template <typename T>
void show_v(vector<T> &v)
{
    for (auto i : v)
    {
        cout << i << " ";
    }
    cout << endl;
}
// 获取列最大最小
void get_min_max(vector<vector<double>> &v, vector<double> &min_, vector<double> &max_, float fac = 10)
{
    // fac 控制界限
    for (auto it : v)
    {
        for (int row = 0; row < 3; row++)
        {
            if (min_[row] > it[row])
            {
                min_[row] = it[row];
            }
            if (max_[row] < it[row])
            {
                max_[row] = it[row];
            }
        }
    }
    for (int i = 0; i < 3; i++)
    {
        min_[i] -= fac;
        max_[i] += fac;
    }
}
// 重载
bool operator<(const std::vector<double> &a, const std::vector<double> &b)
{
    int size = a.size() < b.size() ? a.size() : b.size();
    for (int i = 0; i < size; ++i)
    {
        if (a[i] > b[i])
            return false;
    }
    return true;
}
// 边界内判断
bool in_bord(vector<double> &cor, vector<double> &cmin, vector<double> &cmax)
{
    if (((cmin[0] < cor[0]) && (cor[0] < cmax[0])) &&
        ((cmin[1] < cor[1]) && (cor[1] < cmax[1])) &&
        ((cmin[2] < cor[2]) && (cor[2] < cmax[2])))
        return true;
    else
        return false;
}

void get_ohh(vector<vector<double>> &d,
             vector<vector<double>> &ohh_list,
             vector<vector<double>> &O_ohh)
{
    for (auto itd : d)
    {
        vector<int> idx;
        for (int i = 0; i < ohh_list.size(); i++)
        {
            float dist_ = dist(itd, ohh_list[i]);
            if (dist_ < 10)
            {
                idx.push_back(i);
                O_ohh.push_back(ohh_list[i]);
            }
        }
        // 移除满足距离条件的水分子
        if (!idx.empty())
        {
            for (vector<int>::reverse_iterator it = idx.rbegin(); it != idx.rend(); it++)
            {
                ohh_list.erase(ohh_list.begin() + *it);
            }
        }
    }
}
// 获取公共水分子
void get_com_ohh(vector<vector<double>> &c1_ohh,
                 vector<vector<double>> &c2_ohh,
                 vector<vector<double>> &common_ohh)
{
    for (auto i : c1_ohh)
    {
        for (auto j : c2_ohh)
        {
            bool flag = true;
            for (int k = 0; k < 9; k++)
            {
                if (fabs(i[k] - j[k]) > 1e-6)
                {
                    // 不相等
                    flag = false;
                    break;
                }
            }
            if (flag)
            {
                // 相等
                common_ohh.push_back(i);
            }
        }
    }
}

// 数据计算部分
vector<double> calc(vector<vector<double>> &v)
{
    // 板坐标
    // d1 = data[:248]  // 57.7773 21.1964 39.2873
    // cout << v[0][2] << endl;
    vector<vector<double>> d1(v.begin(), v.begin() + 248);
    // d2 = data[248:496]
    vector<vector<double>> d2(v.begin() + 248, v.begin() + 496);
    // # 两块C板的COO坐标  如果 9个COO 则会有 27行数据  每行为一个原子
    // d1COO = data[221:248]
    vector<vector<double>> d1COO(v.begin() + 221, v.begin() + 248);
    // d2COO = data[469:496]
    vector<vector<double>> d2COO(v.begin() + 469, v.begin() + 496);
    // # H2O 坐标
    // dohh = data[496:]
    vector<vector<double>> dohh(v.begin() + 496, v.end());
    // # 取出 C板的 最大xyz  和最小xyz  增加最大最小值限定  ±10埃米
    // C1_min C1_max C2_min C2_max
    vector<double> C1_min = {9999999.f, 9999999.f, 9999999.f};
    vector<double> C1_max = {-9999999.f, -9999999.f, -9999999.f};
    vector<double> C2_min = {9999999.f, 9999999.f, 9999999.f};
    vector<double> C2_max = {-9999999.f, -9999999.f, -9999999.f};
    get_min_max(d1, C1_min, C1_max);
    get_min_max(d2, C2_min, C2_max);
    // # 取出 C板的COO 最大xyz  和最小xyz  增加最大最小值限定  ±10埃米
    // C1COO_min C1COO_max C2COO_min C2COO_max
    vector<double> C1COO_min = {9999999.f, 9999999.f, 9999999.f};
    vector<double> C1COO_max = {-9999999.f, -9999999.f, -9999999.f};
    vector<double> C2COO_min = {9999999.f, 9999999.f, 9999999.f};
    vector<double> C2COO_max = {-9999999.f, -9999999.f, -9999999.f};
    get_min_max(d1COO, C1COO_min, C1COO_max);
    get_min_max(d2COO, C2COO_min, C2COO_max);
    vector<vector<double>> ohh_list1;
    vector<vector<double>> ohh_list2;
    vector<vector<double>> ohh_list1COO;
    vector<vector<double>> ohh_list2COO;
    vector<double> cor(9, 0);
    for (int i = 0; i < dohh.size() / 3; i++)
    {
        vector<double> cor(9, 0);
        for (int j = 0; j < 3; j++)
        {
            cor[j * 3] = dohh[i * 3 + j][0];
            cor[j * 3 + 1] = dohh[i * 3 + j][1];
            cor[j * 3 + 2] = dohh[i * 3 + j][2];
        }
        // show_v(cor);
        // 剔除距离超出边界的
        // C
        if ((C1_min < cor) && (cor < C1_max))
        {
            ohh_list1.push_back(cor);
        }
        if ((C2_min < cor) && (cor < C2_max))
        {
            ohh_list2.push_back(cor);
        }
        // COO
        if ((C1COO_min < cor) && (cor < C1COO_max))
        {
            ohh_list1COO.push_back(cor);
        }
        if ((C2COO_min < cor) && (cor < C2COO_max))
        {
            ohh_list2COO.push_back(cor);
        }
    }
    // # 获取 C1板附近的±10的水分子
    vector<vector<double>> C1_ohh;
    get_ohh(d1, ohh_list1, C1_ohh);
    // # 获取 C2板附近的±10的水分子
    vector<vector<double>> C2_ohh;
    get_ohh(d2, ohh_list2, C2_ohh);
    // # 获取 C1板COO附近的±10的水分子
    vector<vector<double>> C1COO_ohh;
    get_ohh(d1COO, ohh_list1COO, C1COO_ohh);
    // # 获取 C2板COO附近的±10的水分子
    vector<vector<double>> C2COO_ohh;
    get_ohh(d2COO, ohh_list2COO, C2COO_ohh);
    // # 计算C板 间公共水分子
    vector<vector<double>> C1C2_common_ohh;
    get_com_ohh(C1_ohh, C2_ohh, C1C2_common_ohh);

    // # 计算COO 间公共水分子
    vector<vector<double>> C1C2COO_common_ohh;
    get_com_ohh(C1COO_ohh, C2COO_ohh, C1C2COO_common_ohh);
    // cout << "C1----" << C1_ohh.size() << endl;
    // cout << "C2----" << C2_ohh.size() << endl;
    // cout << "C1COO-" << C1COO_ohh.size() << endl;
    // cout << "C2COO-" << C2COO_ohh.size() << endl;
    
    // common_ohh[0] = C1C2_common_ohh.size();
    // common_ohh[1] = C1C2COO_common_ohh.size();
    // show_v(common_ohh);
    vector<double> common_ohh;
    common_ohh.push_back(C1C2_common_ohh.size());
    common_ohh.push_back(C1C2COO_common_ohh.size());
    mu.lock();
    cout << "C-C-com-ohh : " << C1C2_common_ohh.size() << "---" << C1C2COO_common_ohh.size() << endl;
    mu.unlock();
    // C1板附近满足条件的水分子: 592 个
    // C2板附近满足条件的水分子: 618 个
    // C1-C2板附近公共水分子: 24 个
    // C1板COO附近满足条件的水分子: 551 个
    // C2板COO附近满足条件的水分子: 570 个
    // C1COO-C2COO板附近公共水分子: 23 个
    return common_ohh;
}


// 逐行读取
void readTxt(string file)
{
    ifstream infile;
    infile.open(file.data()); //  将文件流对象与文件连接起来
    assert(infile.is_open()); //  若失败,则输出错误消息,并终止程序运行
    // 所有坐标  总行数减去文件头
    vector<vector<vector<double>>> data;
    for (int it = 0;; it++)
    {
        bool flag = false;
        vector<vector<double>> v(LINES - H_LINE, vector<double>{0.f, 0.f, 0.f});
        string s;
        for (int i = 0; i < LINES; i++)
        {
            if (getline(infile, s))
            {
                if (i >= H_LINE)
                {
                    split(s, v[i - H_LINE], 0);
                }
            }
            else
            {
                flag = true;
                break;
            }
        }
        if (flag)
            break;
        else
            data.push_back(v);
    }
    // show_info(v);
    infile.close(); //关闭文件输入流
    // 获取核心数目
    int core = thread::hardware_concurrency();
    // 记录公共水分子数目
    vector<vector<int>> common_ohh(data.size());
    // 锁

    // 多线程逻辑----没有返回值
    // for (int it = 0; it < (data.size() / core) + 1; it++)
    // {

    //     thread *pth(new thread[core]);
    //     for (int i = 0; i < core; i++)
    //     {
    //         if (it + i >= data.size())
    //             break;
    //         else
    //             pth[i] = thread(calc, ref(data[it + i]));
    //             //pth[i] = thread(calc, ref(data[it + i]), ref(common_ohh[it + i]));
    //     }
    //     for (int i = 0; i < core; i++)
    //     {
    //         if (it + i >= data.size())
    //             break;
    //         else
    //             pth[i].join();
    //     }
    // }
    // 有返回值
    vector<vector<double>> common_ohhs;
    for (int it = 0; it < (data.size() / core) + 1; it++)
    {
        vector<future<vector<double>>> fu;
        fu.reserve(core);

        for (int i = 0; i < core; i++)
        {
            if (it + i >= data.size())
                break;
            else
                fu.push_back(async(calc, ref(data[it + i])));
        }
        for (int i = 0; i < core; i++)
        {
            if (it + i >= data.size())
                break;
            else
                common_ohhs.push_back(fu[i].get());
        }
    }
    show_info(common_ohhs);
}

int main()
{
    clock_t total_time = clock();
    readTxt("1.atom"); // 98901 100695
    cout << (clock() - total_time) / 10 << endl;
}
